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Abstract 

The family of generalized Pseudo-Spectral Time Domain (including the Pseudo- 
Spectral Analytical Time Domain) Particle-in-Cell algorithms offers substantial 
versatility for simulating particle beams and plasmas, and well written codes 
using these algorithms run reasonably fast. When simulating relativistic beams 
and streaming plasmas in multiple dimensions, they are, however, subject to the 
numerical Cherenkov instability. Previous studies have shown that instability 
growth rates can be reduced substantially by modifying slightly the transverse 
fields as seen by the streaming particles . Here, we offer an approach which com¬ 
pletely eliminates the fundamental mode of the numerical Cherenkov instability 
while minimizing the transverse field corrections. The procedure, numerically 
computed residual growth rates (from weaker, higher order instability aliases), 
and comparisons with WARP simulations are presented. In some instances, 
there are no numerical instabilities whatsoever, at least in the linear regime. 

Keywords: Particle-in-cell, Pseudo-Spectral Time-Domain, Relativistic beam, 
Numerical stability. 


1. Introduction 

The Pseudo-Spectral Time Domain (PSTD) Particle-in-Cell (PIC) algorithm 
advances Fourier-transformed electromagnetic fields in time according to the 
difference equations [T], 

E n+1 = E n _ ik x B »+y2 At _ jn+Va Atj ( 1 ) 

B "+ 3 / 2 = B n+1 / 2 + ik x E n+1 Af. (2) 

Simulation particles, of course, require fields and produce currents in real, not 
Fourier, space. Consequently, Fourier transforms of both fields and currents are 
performed at each time step. The fields and currents in real space typically are 
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located at the nodes of a regular multidimensional mesh, with magnetic fields 
B and currents J offset a half time step from electric fields E. 

The PSTD algorithm has been generalized by the authors in two ways mm- 
First, the components of k in Eqs. 0 and 0 are replaced by the Fourier 
transforms of various order finite difference approximations to spatial derivatives 
on a grid [3], which has the effect of narrowing the effective width of the particle 
interpolation stencils. For instance, ki might be replaced by 

sin (kiAxi/2) 

Axi/2 ’ 


yielding a pseudo-spectral realization in k-space of the standard Finite-Difference 
Time Domain (FDTD) algorithm. (In this context the components of k itself 
can be viewed as an infinite order approximation.) Second, the field solver 
is modified to increase the relatively small PSTD Courant time step limit by 
a factor of N, an integer. This modification is equivalent mathematically to 
sub-cycling the field solver N times but is faster computationally. Note that 
the generalized PSTD algorithm reduces to Haber’s Pseudo-Spectral Analytical 
Time-Domain (PSATD) algorithm |5] in the limit of infinite N. Both the details 
of and the motivation for the generalized PSTD algorithm, dubbed the Pseudo- 
Spectral Arbitrary Order Time Domain (PSAOTD) algorithm, are described 
thoroughly in [5]. 

The numerical Cherenkov instability ia 0 is observed commonly in PIC 
simulations of relativistic particle beams and streaming plasmas; e.g., MGS], 
where it can be fast growing and strongly disruptive. At a minimum the instabil¬ 
ity increases beam emittance m- It arises from nonphysical coupling between 
the usual electromagnetic modes, possible distorted by numerical effects, and 
spurious beam modes na da. The dispersion relation describing the numerical 
Cherenkov instability has the general form, 


Co+Wp Y, csc 

m z 

+ 7 _2 Wp (Y c 3z CSC 


(w - k' z v) ^ 


(w - k' z v) 


+U} P Y (C2* + 7 2 C 2 z) esc 
At 
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(w - k' z v) 


At 
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where the C 1 ,; are complicated functions of the instability frequency u>, wave 
numbers k, time step At, cell size, and alias index m z . The particular form 
of the Ci is determined by the details of the numerical algorithm, given in [ 2 j 
for PSAOTD. The relativistic beam itself is characterized by its normalized 
energy 7 , axial velocity v = (l — 7 ^ 2 ) ^, and normalized relativistic plasma 
frequency w p . Of particular importance are the resonances at w = k' z v , with 
alias wave numbers k' z = k z + m z 2n/Az. (The beam propagates along the z- 
axis.) Each alias can trigger an instability, with the most rapid growth typically 
occurring for m z = 0,-1. High order interpolation, say cubic, significantly 
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reduces growth rates associated with higher order aliases. Although peak growth 
rates typically occur at large wave numbers, nontrivial growth rates can occur 
at quite small wave numbers, especially for m z = o mug. The nonresonant 
m z = 0 instability has two branches ES EH, the primary branch associated 
with the C 2 X term of the dispersion relation, and the secondary branch with 
the C 1 term. The former occurs over a wide range of wave numbers, while 
the latter occurs at k z ss 2 Az and small k x , if at all. See Fig. [lj Unless 
otherwise noted, parameters for these and other figures are oj p = 1, 7 = 130, and 
Aa; = A z = 0.3868. The k-space grid is 65x65, corresponding to a simulation 
grid of 128x128. 

Until recently, the numerical Cherenkov instability was ameliorated by a 
combination of digital filtering, numerical damping, higher order interpolation 
(typically cubic), and astute parameter choices [S]|T5j. Papers within the past 
few years have offered alternative options, often entailing minor corrections to 
the transverse electric and magnetic fields as seen by the particles MEJE]. 
Specifically, transverse electric and magnetic fields as they are interpolated to 
the particles are multiplied by k-dependent correction factors and *Fb that 
vary as 1 + O (fc 2 ) for small k. The distinction between such correction factors 
and the usual digital filters is made clear in mm- The particular analytical 
correction factors presented in these references are quite effective at ameliorating 
the first branch of the m z = 0 instability but ineffective at ameliorating the 
second, for reasons presented later in this article. 

As an alternative, this paper numerically computes correction factors *F e 
and ’Fb just sufficient to completely eliminate both branches of the m z = 0 
numerical Cherenkov instability. The process is straightforward in principle. 
The m z = 0 dispersion function is akin to a fifth order polynomial with real 
coefficients. If it crosses the w-axis in five places (i.e., five real roots), all roots 
are stable. However, if parameters are varied such that the curve crosses the 
axis in only three places or one, then the dispersion function has two or four 
complex roots, respectively, the remaining roots being real. The transition 
occurs where the curve becomes tangential to the w-axis, i.e., where both the 
dispersion function and its derivative with respect to w vanish. This determines 
5 'e and 'Fb. In fact, it can be shown that the m z = 0 portion of Eq. ^ can 
be represented as a fifth order polynomial, although this is not required for the 
argument of this paragraph to be true. 

Fig. i computed for vAt /Az = 0.9 and k z = k x = 5.20, is illustrative. The 
left plot displays the dispersion function, Eq. ©> with m z = 0 only, for the 
complete frequency range, —^/a t < w < At. Shown are the results for the 
uncorrected case, >Fb = ’F b — 1 (labeled Base), for the C^x correction factors 
described in [2] (labeled C- 2 x), and for the numerically determined optimal cor¬ 
rection factors (labeled Opt). (The C^ x correction factors are chosen so that the 
term C^ x vanishes at w = k z v). The three curves are essentially indistinguish¬ 
able on this scale, with one crossing at the far left and two or four crossings 
at the right. The right plot displays a blow-up of the critical frequency range. 
Each curve has one crossing in this range, and a second crossing is off-scale to 
the right. However, the Opt curve is seen also to be tangent to the w-axis at a 


3 


point. Hence, it is stable, and the other two cases are not. The corresponding 
m z = 0 instability growth rates are 0.0749, 0.0014, and 0 for Base, C 2 X , and 
Opt. 

The next section provides sample numerical correction factors and the pro¬ 
cedure used to obtain them. The third section then provides corresponding nu¬ 
merical growth rates from Eq. (|3j, along with corroborating instability growth 
rates from WARP two-dimensional simulations m- a short concluding section 
completes the paper. 

The analytical and numerical linear growth rate analyses were performed 
using Mathematica m- 

2 . Numerical Procedure for Obtaining \Fe and \Fb 

The factors 4 /e and 4 /b appear linearly in the m z = 0 dispersion function, 
obtained from Eq. [3] by omitting all but the m z = 0 terms in the summations, 
and in its derivative with respect to w. As a consequence, the dispersion function 
and its derivative together can be inverted readily to obtain ’I'e and as 
functions of the tangent point, to, although the actual expressions are quite 
complicated. 

Fig. ! (left) is a typical parametric plot of $5 vs. 't® as u is varied. Its 
parameters are identical to those of Fig. [2j The region below the diagonal 
curve is stable, above it unstable. (Despite appearances, it bends slightly and 
does not pass through the origin.) Although any choice of 4/^ and along 
or below this curve stabilizes the m z = 0 mode, a pair close to ( 1 , 1 ) is to 
be preferred for minimizing numerical impact on the physics to be simulated. 
Further constraining the choice by I'b < 1, < 1 suggests the following 

procedure: 

1. If numerical instability growth is zero at (1,1), set ’Fg = 1, ’F# = 1. 

2. Otherwise, search the lines b = 1 and \Fe = 1 for the point nearest (1,1) 
at which the dispersion function and its derivative simultaneously vanish. 
Usually, that point will be on the 4/^ = 1 line, as illustrated in the left 
plot of Fig. [3] 

3. Otherwise, search for the point nearest (1,1) at which the dispersion 
function and its derivative simultaneously vanish in the interior of the 
0<4's<l A 0<4'£’<1 region. It is fortunate that this situation is 
rare, because such points are expensive to find computationally and easy 
to miss. 

4. Otherwise, set 4^ = ^e = 0. The numerical instability cannot be elimi¬ 
nated, short of excising that portion of k-space. 

For completeness, Fig. [ 3 ] (right) provides the parametric plot of 4 /b vs. 4>e for 
k z = 4.697, k x = 0.635, which lies at the center of the second branch of the 
m z = 0 numerical Cherenkov instability in Fig. [T] Stable regions lie above the 
upper diagonal curve and below the lower diagonal curve. Only a very small 
reduction from unity in e is sufficient to eliminate the numerical instability. 
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Plots of 4 <e and V F a for the parameters of Fig. [2] are shown in Fig. [4] with 
the m z = 0 resonance curve superimposed. *P e is equal to unity everywhere but 
(1) in the small region, barely visible at the bottom center of the plot, where the 
second branch of the m z = 0 numerical Cherenkov instability appears in Fig. 
[lj (2) near the resonance curve, where 4/# = 0, and (3) at three isolated points 
where = 0 also, 4^ — 0.97 for the second branch of the m z = 0 numerical 
instability. In contrast, < 1 everywhere except those k values where the 
m z = 0 growth rate is zero anyway. Note that 1 — 4^ ~ 0 over a wide range of 
small k. 

The three isolated = 4^ = 0 points appearing in both plots are asso¬ 
ciated with a weak instability predicted by Birdsall and Langdon in Problem 
5-9a of their well known text m■ Because this instability has a very narrow 
bandwidth in k-space, it appears only where the instability band happens to 
overlap nodes on the k-space numerical grid. Analysis of this Courant-like in¬ 
stability for the generalized PSTD algorithm is provided in [2], and analysis plus 
growth rates for the PSATD algorithm in m ■ This instability appears to be 
unimportant in practice except, perhaps, at low 7. 

Applying the correction factors of Fig. [4] for the parameters used to obtain 
Fig. 0 completely eliminates m z = 0 numerical instability growth but has 
minimal effect on higher order aliases except for k-space regions where 4'g or 
4 /E are much less than unity. Hence, we choose to truncate fields in k-space 
according to 


k > a min 


7T 7r 

A z ’ v At 


(4) 


in order to minimized peak growth rates of the first branch of the m z = — 1 
instability. Fig. [5] (left) displays the resulting growth rates for a = 0.85. Trun¬ 
cating k-space according to Eq. Q reduces the peak growth rate from 0.168 
to 0.071, and the Fig. [4] correction factors further reduce it to 0.026. Most of 
the residual growth is associated with the second m z = — 1 branch. The second 
branch of the m z = 0 numerical Cherenkov instability is absent, as expected. In 
contrast, Eq. @ with C 2 X correction factors reduces peak growth comparably, 
but with nontrivial growth rates spread over a larger range in k-space; see Fig. 
[5] (right). In particular, it causes the second branch of the m z = 0 numerical in¬ 
stability to shift toward k x = 0 but does nothing to eliminate it. This is because 
the topology of the vs. 4 /e parametric plot changes rapidly for very small 
k x near this second instability branch. As a consequence, slightly reducing >F e 
there actually can create an instability where none otherwise exists. 


3. Numerical Instability Growth Rates 

Although k-space truncation with a = 0.85 yields reasonably effective in¬ 
stability suppression, a = 0.60 is much better. Additionally, choosing a = 0.60 
better accommodates comparisons with earlier papers DU EG]. Fig. [6] (left) 
depicts peak growth rates for N = 1,2, 3, 4, 8, 00, the last being PSATD. Peak 
growth rates are relatively insensitive to N. As always, the rapid variation of 
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peak growth with v At /az reflects the rapidly changing locations of higher order 
resonances on the k-space grid. Note that the largest growth rates are associ¬ 
ated with ” Ai /Az 0.9 , used as the sample case in the preceding Section. No 
growth whatsoever is predicted for ' At /Az = 1.0, because all aliases coincide 
there. Peak growth rates for N = 4 and various order finite difference approxi¬ 
mations to k are shown in Fig. [ 6 ] (right). Results are similar to those in the left 
plot. Superimposed on the growth rate curves are results from several WARP 
simulations m- Agreement is very good. In addition, a simulation was run 
to t = 10000 for ” Ai /az = 1.0 with no sign of instability growth. However, at 
still larger values of time, instability growth at a rate of about 0.003 developed 
at wave numbers consistent with the second branch of the m z = 0 numerical 
Cherenkov instability. Perhaps, a slight change in the beam distribution func¬ 
tion over time moved this instability above threshold. If so, reducing T e very 
slightly there might prevent this very late time growth. 

The procedure introduced in this paper remains robust at low 7 . Fig. [7] 
presents peak growth rates for 7 = 3, a = 0.8, N = 2, 3, 4, 8 , 00 , and infinite 
order. Results are comparable to the best obtained earlier; see Fig. 7 of EH. 
Both here and in m, E z is offset by one-half cell to suppress the usual quasi¬ 
electrostatic numerical instability that occurs at low 7 |21[ [22| 123] . The peak in 
Fig. 0 at ,;At /Az = 0.9 is larger than that in Fig. [ 6 ] only because a is larger. 

4. Conclusion 

New field correction factors have been derived that completely eliminate the 
m z = 0 numerical Cherenkov instability for the generalized PSTD algorithm, 
including the PSATD algorithm. When combined with a sharp cutoff digital 
filter at large k, these correction factors reduce peak growth rates to less (often 
much less) than 0.01 of the beam’s relativistic plasma frequency. These coef¬ 
ficients are optimal in the sense that they differ from unity only slightly over 
a broad portion of k-space while eliminating both m z = 0 numerical instabil¬ 
ity branches. A disadvantage from the implementation perspective is that the 
coefficients must be computed numerically for each set of simulation param¬ 
eters. Software for doing so is available in Mathematica CDF format [2f|| at 
http: //hifweb.lbl.gov/public/BLAST / Godfrey/ 
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Figure 1: Growth rates and resonance curves from PSAOTD dispersion relation for m z = 
[—2, +2], vAt/Az = 0.9, and N = 4. The secondary m z = 0 branch, although slower growing 
and confined to a small region in /c-space, is inconveniently located at small k x . 
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Figure 2: Dispersion function for Base, C^x, and Opt correction factors. Left, entire range of 
frequencies (curves indistinguishable at this scale). Right, narrow range of frequencies near 
second crossing point. Parameters are as in Fig. [T] with k z = k x = 5.20. 
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Figure 3: *$/b — ^E parametric plots for the parameters of Fig. [T] left, k z = 5.20, k x = 5.20, 
and right, k z = 4.696, k x = 0.635. Shaded regions are stable. 
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Figure 4: 4/# (left) and 4 /b (right) correction terms for the parameters of Fig. [T] Superim¬ 
posed at the far upper right of each plot is the m z = 0 resonance curve. 
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Figure 5: Growth rates for Fig. 0 parameters, optimal (left) or Cix (right) correction factors, 
and a = 0.85. 
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Figure 6: Peak numerical instability growth rates vs. vAt//± z with for various N (left) and 
for various order finite difference approximations to k (right); a = 0.6. 
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Figure 7: Peak numerical instability growth rates vs. ' Ai/aj with 7 
a = 0 . 8 . 
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